#include "../header.h"

#define alpha 1
#define beta 1.05
#define gamma 1.03
#define c 1.04
#define Dn 0.005
#define Dm 0.005

double RabbitStartDensity(double x,
                          double y) {
  if (x <= a || x >= b || y <= a || y >= b)
    return 0;
  else
    return 100;
}

double RabbitsInSquare(double x,
                       double y) {
	double temp;
	if(x == 0 || x == 1 || y == 0 || y == 1) return 0;
	temp = -0.08*((x-0.15)*(x-0.15)+(y-0.15)*(y-0.15))+1;
	if(temp < 0) return 0;
		else return temp;
}

double WolfStartDensity(double x,
                        double y) {
  if (x <= a || x >= b || y <= a || y >= b)
    return 0;
  else
    return 100;
}

double WolfsInSquare(double x,
                     double y) {
  double temp;
	if(x == 0 || x == 1 || y == 0 || y == 1) return 0;
	temp = -0.1*((x-0.15)*(x-0.15)+(y-0.15)*(y-0.15))+1.25;
	if(temp < 0) return 0;
		else return temp;
}

void GetInitialDistribution(double **M, double **N,
                            int from_point, int to_point) {
  int i,j,points,ypoints;

  ypoints = ysteps;
  points = to_point-from_point+1;  
  for(i=0;i<ypoints;++i)
    for(j=0;j<points;++j) {
      N[i][j] = RabbitsInSquare((from_point+j)*hx,i*hy);
      M[i][j] = WolfsInSquare((from_point+j)*hx,i*hy);
    }
}

//returns d2f/dx2
double SecondDif(double **x, int m, int n, char DifVar) {
//  int mmax,nmax;
//  mmax = (int) ((double)(b-a))/hy;
//  nmax = (int) ((double)(b-a))/hx;
  switch (DifVar) {
    case 'x':
      return (x[m][n+1]-2*x[m][n]+x[m][n-1])/((double)hx*hx);
      break;
    case 'y':
      if ((m == 0) || (m == ysteps - 1)) 
        return x[m][n];
      else 
        return (x[m+1][n]-2*x[m][n]+x[m-1][n])/((double)hx*hx);
      break;
  }
}

/*
  Below we solve system:
              dN(t)/dt = (alpha - c*M(t))*N(t) == f1(M,N)
              dM(t)/dt = (-beta + gamma*N(t))*M(t) == f2(M,N)
*/

double Getf1(double M, double N) {
  return (double) (alpha - c*M)*N;
}

double Getf2(double M, double N) {
  return (double) (-beta + gamma*N)*M;
}

void GetNextNonLinearIteration(double *Mnext, double *Nnext, double Mprev, double Nprev) {
  double k1N = tay*Getf1(Mprev,Nprev);
  double k1M = tay*Getf2(Mprev,Nprev);
  double k2N = tay*Getf1(Mprev+k1N/2,Nprev+k1M/2);
  double k2M = tay*Getf2(Mprev+k1N/2,Nprev+k1M/2);
  double k3N = tay*Getf1(Mprev+k2N/2,Nprev+k2M/2);
  double k3M = tay*Getf2(Mprev+k2N/2,Nprev+k2M/2);
  double k4N = tay*Getf1(Mprev+k3N,Nprev+k3M);
  double k4M = tay*Getf2(Mprev+k3N,Nprev+k3M);

  *Mnext = Mprev + (k1M+k2M+k3M+k4M)/6;
  *Nnext = Nprev + (k1N+k2N+k3N+k4N)/6;
}

/* 
  And now the system is
  (Ntemp-Nprev)/tay = (alpha - c*M(t))*N(t)
  (Nnext-Ntemp)/tay =  Dn*laplassian(N)
  (Mtemp-Mprev)/tay = (-beta + gamma*N(t))*M(t)
  (Mnext-Mtemp)/tay =  Dm*laplassian(M)
*/

void GetNextIteration(double **MPrevGrid, double **NPrevGrid,
                      double **MNextGrid, double **NNextGrid,
                      int from_point, int to_point) {
  double Ntemp,Mtemp;
  int i,j,points;

  points = to_point-from_point+1;
  for(i=0;i<(int)ysteps;++i)
    for(j=1;j<=(int)points;++j) {
      Ntemp = 0;
      Mtemp = 0;
      GetNextNonLinearIteration(&Mtemp,&Ntemp,MPrevGrid[i][j],NPrevGrid[i][j]);
      MNextGrid[i][j] = Mtemp + (double)tay*Dm*(SecondDif(MPrevGrid,i,j,'x')+SecondDif(MPrevGrid,i,j,'y'));
      NNextGrid[i][j] = Ntemp + (double)tay*Dn*(SecondDif(NPrevGrid,i,j,'x')+SecondDif(NPrevGrid,i,j,'y'));
    }
}
/*
void GetNextIteration(double **MPrevGrid, double **NPrevGrid,
                      double **MNextGrid, double **NNextGrid,
                      int from_point, int to_point) {
  int i,j,points;

  points = to_point-from_point+1; 
  for(i=0;i<ysteps;++i)
    for(j=1;j<points;++j) {
      MNextGrid[i][j] = MPrevGrid[i][j];
      NNextGrid[i][j] = NPrevGrid[i][j];
    }
}
*/      












